In this practical, we will analyse data pertaining to customer waiting times in a governmental office. The data provided to us is the \(waiting.csv\) file are daily average waiting times from January 3, 2017 to October 26, 2018.

Our goal is to model and monitor the time customers spend waiting, and to take action if the service quality deteriorates too badly, i.e. if waiting times are unacceptably high.

Exercise A

Read in the data and plot the daily waiting times. Compute and display the five-number summary (boxplot) of the daily waiting times for each weekday. What do you observe?

Daily plot

Boxplot

Table 1: Five-number summary by weekday
Weekday Min Q1 Mean Q3 Max
Friday 148 342.25 485 622.00 918
Monday 116 297.00 446 558.75 999
Thursday 81 217.00 373 487.25 938
Tuesday 90 262.00 431 550.00 1154
Wednesday 95 293.25 471 597.00 1182

\(Conclusion\) : our main observations regarding the charts below are the following:

  • The daily plot chart show a continue oscillation of the data during the 447 observations (250 obs in 2017 and 197 in 2018).
  • There are high deviation values regarding the mean for each day.
  • The average waiting time per weekday are similar within them, nevertheless, the highest average is made on friday.
  • The highest waiting time is on wednesday.
  • The extreme values are located in the first 4 weekdays and in general they are higher from 900 seconds.

Exercise B

Using a normal approximation produce an upper control limit or confidence line for the waiting times at the one-year return level. What do you notice? Using a Q-Q plot, comment on whether the Normal model is appropriate.

Upper Control Limit

########Assuming that waiting times are normally distributed

#95% CI--> alpha = 0.05 We can get z(alpha/2) = z(0.025) from R:
nrow<- nrow(data)
qnorm<- qnorm(.975)   #1.959964
meanData<- mean(data$Average.Wait.Seconds) #440.6398
sdData<- sd(data$Average.Wait.Seconds) #206.19

###########Our margin of error is 

upper_control_limit<- mean(data$Average.Wait.Seconds)+
                            qnorm*(sdData/sqrt(nrow))

Probability of exceeding

#If X~N(μ, σ)--> pnorm(x, μ, σ) function to calculate P(X > upperbound).
# r=1

probability_exceed_UCL<- 1 - pnorm( upper_control_limit, 
                         meanData, sdData)
probability_exceed_UCL
## [1] 0.4630697

QQ-Plot

\(Conclusion\): our main observations regarding the analisis are:

  • The upper bound with a normal approximation is very close the mean of the waiting time (Upper bound= 460, Mean= 441) but with this analysis we do not consider information on the behavior of extreme values, which is key to reducing risk –> central limit theorem.
  • 46% of the time we will exceed the upper limit.
  • In the QQ-plot we can confirm that the data is not normal distributed in the limits, because it is right-skewed (or positively skewed) due to the desviation of the straight line in both extrems.

Exercise C

The aim of this analysis is to focus on negative effects (long wait times). Explain how you would proceed using 1) a block maxima and 2) a peaks-over-threshold approach. For each method, carry out the required data aggregation and transformation.

Block Maxima approach

We have proposed the following steps:

  • First, we need to obtain our \(M_{n,i}\), where \(n\) is the length of the block–> in this case we have choosen weekly because we only have 2 years of data, that is quite short for analyzed by year.

Block Maxima: \[M_{i}=max\{X_{(i−1)n+1},\cdots,X_{in}\},i=1,⋯,m\]

  • Then, we fit the a Generalized extrem Value, assuming that it follows \(GEV(\mu_{n},\sigma_{n},\xi)\) distribution and we will focus on maximize the log-likelihood and the estimation of the parameters.

  • In this step, we decide to test if the parameters correspond to the Gumbel distribution.

  • Finally, we asses the model and compute m-year return levels.

Divide data

#we should first extract the maximum temperature values-->  block maxima.

#week variable creation
# derive AMS for maximum precipitation

data_ts <- xts::xts(data$Average.Wait.Seconds, order.by = data$Date) #time series creation
wms_data <- xts::apply.weekly(data_ts, max) #create the weekly maximal serie

Plot

Fit to GEV

#According to the Fisher–Tippett–Gnedenko theorem, the distribution of block maxima can be approximated by a generalized extreme value distribution.

gev.fit2 <- extRemes::fevd(as.vector(wms_data),method = "MLE", type = "GEV")

Testing Parameters

(#tab:table_parameters )This table show the parameters of the GEV fit model
Location Scale Shape
Estimates 585.94 170.64 -0.13
Std.errors 19.87 14.06 0.07 
######Testing the Gumbel distribution (ξ=0) hypothesis with LRT#####################
#################################################################################

gev_gumb <- fevd(as.vector(wms_data), type="Gumbel", units="deg C")

test_gumbel<- lr.test(gev_gumb,gev.fit2)

test_gumbel
## 
##  Likelihood-ratio Test
## 
## data:  as.vector(wms_data)as.vector(wms_data)
## Likelihood-ratio = 2.478, chi-square critical value = 3.8415, alpha =
## 0.0500, Degrees of Freedom = 1.0000, p-value = 0.1154
## alternative hypothesis: greater
#leads us to not reject the Gumbel hypothesis !
########################################################################

Plot

\(Conclusion\) : our main observations are the following:

  • Even if the shape value is negative it is small enough to fit a Gumbel distribution as we test with the measure of likelihood-ratio test.

we are going to perform the analysis of the plot of the fit in the question D.

Peaks-over-theshold approach

We have proposed the following steps:

  • First, we need stablish a large enough treshold \(u\), where the following formule need to be fullfill: \[\mathbb{P}\left(X_{i}-u>y|X_{i}>u\right)\] In order to do that, we are going to run a mean residual plot and then we will select the value where the plot start looks linear.

  • Then, we fit the all the dataset to a Generaliaze Pareto distribution.

  • Finally, we asses the model and compute m-year return levels.

Setting the threshold

Plot data

Fit GP

# maximum likelihood estimation
pot.fit <- fevd(as.vector(data_ts), method = "MLE", 
                type="GP", 
                threshold = threshold)

Plot GP

\(Conclusion\) : our main observations are the following.

  • We have changed the threshold level a couple of times and we chose the value 700 since in the density graph we can see what the model captures better.

we are going to perform the analysis of the plot of the fit in the question D.

Exercise D

Propose a model for the data you have processed in the previous question. Make sure to justify your model choice using residuals and goodness-of-fit tests. (Hint: you may use the evd package.)

Now, we will use likelihood-ratio test for two nested extreme value distribution models,

the function used comes from the extRemes package.

Quantile plot

par(mfrow=c(1,2))
plot(gev.fit2, type = "qq2") 
plot(pot.fit, type = "qq2")

Figure 1: Quantile plot GEV-GP

Quantile plot GEV-GP

\(Conclusion\): our main observations regarding the analisis are the following.

  • As you can see in the charts, the GEV model has a better fit at the beginning, however, the GP model is better in the extreme values.
  • In addittion, the AICs of the GP model is the one that show a better fit of the data. AIC(GEV)= 1242.5 and AIC(GP)= 586.
  • Finally, we want to avoid overfitting that is more common with the first model.

Exercise E

Finally, use your model to derive an upper control limit or confidence line for the waiting times at the one-year return level.

Computing CI

(#tab:selected CI)Confidence interval for one-year return level
name values
one year return level 1140.82
Lower bound 1041.54
Upper bound 1240.09

Plots Bounds

\(Conclusion\) :

  • Meaning that once a year we are going to exceed the value of 1140.82 seconds for the waiting time with a 95% of confidence.

In this second practical, we will look at product sales data belonging to a large American retailer.

As the person in charge of re-stocking products after they are sold, you would like to avoid a ‘stock-out’, that is, when on a given day there are not enough products to satisfy customer demand.

Consider the observations recorded in sales 1.txt, which contain the daily sales volumes of Product 1 over n = 1913 days.

Exercise a

Construct a figure which shows the daily sales volume, and study the values just before the product is out of stock. What do you notice? Explain why this is related to stock-out.

plot(sales$Sales, type = "l", main = "Daily Sales")
outofstock <- which(diff(sales$Sales) == -sales$Sales & sales$Sales !=0)
## Warning in diff(sales$Sales) == -sales$Sales: longer object length is not a
## multiple of shorter object length
abline(v = outofstock, col = "red")

sales %>% dplyr::mutate(date = as.Date(1:1913)) %>% dplyr::filter(Sales == 0)
##   Sales       date
## 1     0 1970-11-28
## 2     0 1971-01-20
## 3     0 1971-01-21
## 4     0 1971-11-29
## 5     0 1973-06-29
## 6     0 1973-06-30
## 7     0 1973-11-28
## 8     0 1974-11-28

It seems there is some seasonality element in the data, as the product is out of stock alomsto every year at the same time (around the end of November for almost every year), hence it could be avoided as it is repeated thorugh different years.

Exercise b

Assume Product 1 is a perishable good sold at a constant price of 10.– over the time period considered. Items of Product 1 are produced at cost 1.– and items which are unsold can be salvaged at 10% of cost, that is, at value s = 10 ct. Suppose that the sales of Product 1 are distributed according to F. Using the newsvendor model, give a formula for the critical fractile q = F−1(p).

Otherwise, we can calculate it by hand.

#price
pr <- 10
#cost of production
c <- 1 
#salvage cost
s <- 0.1
#cost of shortage 
cs <- pr-c
#cost of overage 
co <- c-s
#critical fractile
p <- cs / (cs + co)

In this case, the formula to calculate the critical fractile is the following:

\[ q = F^-1 (p)\]

is equal to

\[ q = F^{-1} (\frac {cs}{cs + co}) \]

The cs as the cost of shortage, which corresponds to the cost per unit it would incur if there is still demand but we are out of stock, hence a missed sale which correspond to the price minus the costs (p - c). On the other hand co is the cost of overage, which is the cost per unit of having too many products compared to the demand, hence a product produced and unsold which is the cost of production minus the salvage cost (c - s).

Exercise c

Explain how the critical fractile is related to the Value at Risk, and give an interpretation of the Expected Shortfall at the level p.

The critical fractile, balances the cost of being understocked and the total costs of being either overstocked or understocked, and it gives the optimal quantity that should be produced to minimizes the costs. It’s the quantile which optimize the probability that the demand is lower than the quantity Q that has been ordered.

The Value at Risk (VaR) denotes a quantile of the distribution, which is the smallest number such that the random variable X that we are studying exceeds a value x with probability lower than 1 - alpha. In other words, with probability alpha, the variable X will be smaller than VaR. Hence, we are looking at the heavy tails of the distribution of the variable.

The interest is when the newsvendor model gives a high quantile to look at (e.g. above 0.9, like in our case in the previous question). What we want to avoid is to go beyond the critical fractile, because in that case we will not choose the optimal quantity, and we will increase the costs and hence decrease the profit, even though this will happen with a low probability (of 1 - q). This is linkable to the Value at Risk as it is also describing an extreme situation in which the variable takes a value that is higher than a certain threshold, and, as we are in the risk analysis field we are studying the losses on investments, we want to avoid these extreme cases, as they mean a high loss.

The Expected Shortfall at level alpha is the average value of X when it is higher than VaR. It is related to VaR by

\[ ES_\alpha (x) = \frac{1}{1 - \alpha} \int_\alpha^{1}q_u (F_X)du = \frac{1}{1 - \alpha} \int_\alpha^{1} VaR_\beta(L)d\beta \]

The Expected Shortfall at level p will then be the value of the sales if they are higher then the VaR, so for the level p, the ES will be the average value of the sales that are higher than the level p that is calculated as the quantile of the critical fractile q.

Exercise d

Fit a Poisson model to the sales volume data by assuming that Yi ∼ Poisson(μ) and using maximum likelihood estimation for μ. (Hint: you may use MASS::fitdistr). According to this model, what is the estimated critical fractile?

poisson <- fitdistr(sales$Sales, "Poisson")
#lambda of the poisson distribution
lambda <- poisson$estimate
#quantile 
q <- qpois(p, lambda)

In the case of a poisson distribution, we find a critical fractile by using the following formula:

\[ q = F^{-1} (p) = \sum_{i = 0}^p \frac {e{-\lambda}\lambda^i}{i!}\]

And it corresponds to 114.

Exercise e

Suppose we model the sales as an extreme value distribution using a peaks-over-thresholds method. Give the Value at Risk at level p under this model.

#Mean Residual Life Plot
extRemes::mrlplot(sales$Sales, main="Mean Residual Life Plot")

#Finding and appropriate threshold 
extRemes::threshrange.plot(sales$Sales, r = c(100, 200), nint = 120)

#let's fix the threshold at 140 as after it the distribution seems to be flat in the second graph and there is the elbow in the first one
u <- 190

plot(sales$Sales, type = "l", main = "Daily Sales") +
abline(h = u, col = "red")

## integer(0)
#Weights for the observations above the threshold chosen
W <-  sales %>% dplyr::filter(Sales > u) %>% dplyr::mutate(W = Sales - u) %>% dplyr::select(W)

#fitting a GPD, method amle as we want the shape parameter positive
gpd <- gPdtest::gpd.fit(W$W, method = "amle")

#extracting the parameters 
#shape
beta <- gpd[1,]
#scale
xi <- gpd[2,]

#setting the alpha for the VaR to the level we have calculated 
alpha <- p

#calculating Fbar with the approximation 
Fbar <- nrow(W) / nrow(sales)

#calculating the values at risk 
VaR <- (u + (beta/xi) * ( ( ( (1 - alpha) / (Fbar) ) ^ (-xi) ) - 1) )

We wanted to use the first two graphs to determine the therhsold that we want to use. What we can see in the first one is that between the value 100 and 150 the curve seems to flatten, which is confirmed by the second graphs, as the parameters are becoming constant already from before 100. However, choosing a threshold between these two values would give too many observations above it, in our opinion, hence we decide to move it a little bit further up. We decided to be safe and to choose a threshold of 190, so that we have a quantile that is around the 80%.

We move on by calculating the weights of the observations that are above this threshold (that are the points in the graph above the red line).

Then, using these observations we fit a GPD and we estimate the parameters, that we will insert in the formula to calculate the VaR.

The resulting value is 189.9589851.

Exercise f

Perform a binomial back-test over the last 300 days in the dataset for the models in d) and e), using a window size of 365 days. What are your conclusions?

To backtest: 1. Split the data and create a test set including the last 300 observations of the dataset 2. Calculate the VaR for each one of the observations inside the test set using a window of 365 days 3. Compare the estimated VaR to the real values and see if there is any violation, meaning that the real value is higher than the VaR for a given day 4. Consider that the violations follow a binomial distribution(300, 1 - p), test the hypothesis that the observed proportion of violations over the test set is different from the expected proportion given by the binomial distribution with a binomial test

Poisson

1. Test set

We will split the data in two subgroups, the first one being the training set which we will use to estimate the VaR will comprehend the observations up to the last 300 days of the data, while the last 300 observations will be the test set to which we will compare the VaR we will find.

salesTest <- sales$Sales[(nrow(sales)-299):(nrow(sales))]

2. Calculate VaR

Then, we calculate the Value at Risk by fitting a poisson distribution on the window preceeding the observation in which we are interested.

We start by creating a function that works in the following way: 1. use the same probability as before (hence 0.91) 2. fit the poisson on the window 3. determines the lambda of the poisson distribution 4. caclulate the VaR by taking the 0.91 quantile of the poisson we fit

poissVaR <- function(df){ 
  
  #definition of the probability 
  p <- 0.91
  
  #fitting the poisson
  poisson <- fitdistr(df$Sales, "Poisson")
  
  #lambda of the poisson distribution
  lambda <- poisson$estimate
  
  #quantile, which corresponds to the VaR
  VaRpoiss <- qpois(p, lambda)
}

Then, we can create a loop to get the VaR for the observations inside the test set.

#create a vector to store the values 
poissVAR <- as.numeric(1:300)

for(i in 1:300){
  #select window 
  wind <- sales %>% dplyr::slice((1913-i):(1913-i-365))
  #calculate VaR
  poissVAR[i] <- poissVaR(wind)
}

poissVAR
##   [1] 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106
##  [19] 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106
##  [37] 106 106 106 106 106 106 106 106 106 106 106 106 106 106 105 105 105 105
##  [55] 105 106 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105
##  [73] 105 105 105 105 105 105 105 105 105 105 105 106 105 105 105 105 105 106
##  [91] 106 106 106 105 105 105 106 106 106 105 105 105 105 106 106 106 105 105
## [109] 105 105 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106
## [127] 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 105 105
## [145] 106 106 106 106 106 106 105 105 105 106 105 105 105 105 105 105 106 106
## [163] 106 105 105 105 105 106 106 105 105 105 105 105 105 105 105 105 105 105
## [181] 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105
## [199] 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105
## [217] 105 105 105 105 105 105 105 105 105 105 105 106 106 106 106 106 106 106
## [235] 106 106 106 107 107 106 107 107 107 107 107 107 106 106 106 106 106 106
## [253] 106 106 106 106 106 106 107 107 107 107 107 107 107 107 107 107 107 107
## [271] 107 107 107 107 107 107 107 107 107 107 107 107 107 107 107 107 107 107
## [289] 107 107 107 107 107 107 107 107 107 107 107 107

As we can see they take they almost take the same values, this is due to the fact that the window is quite large (365), hence the change from one window to the subsequent one will not be very big, as we just drop the furthest away observation to gain the most recent one that is the day before the one we are considering. Hence, the fitting of the poisson distribution will be almost the same, as it will be the lambda and the quantile.

3. Violations

Now, we compare the VaR we estimated to the real values, to see how many violations there are. The should be following a bernoulli distribution with n = 300 and p = 1 - alpha.

#number of observed violations
poiss.violations <- sum(ifelse(poissVAR > salesTest, 0, 1))
#number of theoretical violations, given by the mean of the bernoully distribution
ev <- 300*(1 - p)

We find out that there are 89 violations.

We can compare this value to the theoretical one that we expect to get, which is 27.2727273.

4. Binomial test

Binomial test

We can test whether they are statistically significantly different by using a binomial test, which has the following characteristics:

\[H_0: \pi = \pi_0\] In our case the theoretical π will equal to 0.0909091, while the probability found in the observations is 0.2966667. We want to test if waht it observed is different from the expected value.

(bintest.poiss <- binom.test(poiss.violations, 300, p = 1-p))
## 
##  Exact binomial test
## 
## data:  poiss.violations and 300
## number of successes = 89, number of trials = 300, p-value < 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.09090909
## 95 percent confidence interval:
##  0.2455466 0.3518550
## sample estimates:
## probability of success 
##              0.2966667

We can see that for the poisson, the p-value is really low 4.138435310^{-24}. This means that we will reject the null hypothesis and hence that the distribution of the violations is not following a bernoulli distributions. There has been an underestimation of the quantile we used to calculate the VaR and the values at risk, as we had quite a lot violations compared to the theoretical one. We believe that by increasing the probability in the poisson distribution we would get to a more optimal result.

Peaks Over Threshold

1. Test set

We will split the data in two subgroups, the first one being the training set which we will use to estimate the VaR will comprehend the observations up to the last 300 days of the data, while the last 300 observations will be the test set to which we will compare the VaR we will find.

salesTest <- sales$Sales[(nrow(sales)-299):(nrow(sales))]

2. Calculate VaR

Then we will calculate the VaR for the test set using a moving window of 365 days.

In order to do so, we start by creating a function to calculate the VaR for a given window, which works in the following way: - here we use the 73% quantile to determine the threshold as a higher one would mean too few observations above the line to be able to fit a gpd and it would create too few violations, while a lower one would create the opposite problems, - window of 365 days before the observations that we want to estimate, - calculate the weights for the observations in this window which are above the threshold by taking the differences between the sales and the thershold, - fit a GPD on the weights and estimate the parameters for the shape and the scale, - same alpha as before, which is the p value we calculated, - estimation of Fbar ny dividing the number of observations above the threshold in the window and the total number of observation in the window, - put all the estimators in the formula to calculate the value at risk.

VaRest <- function(df){
  
  #threshold
  u <- 0.73 * max(df$Sales)
  
  #Weights for the observations above the threshold chosen
  W <-  df %>% dplyr::filter(Sales > u) %>% dplyr::mutate(W = Sales - u) %>% dplyr::select(W)
  
  #fitting a GPD, method amle as we want the shape parameter positive
  gpd <- gPdtest::gpd.fit(W$W, method = "amle")
  
  #extracting the parameters 
  #shape
  beta <- gpd[1,]
  #scale
  xi <- gpd[2,]
  
  #setting the alpha for the VaR to the level we have calculated
  alpha <- p
  
  #calculating Fbar with the approximation 
  Fbar <- nrow(W) / 300
  
  #calculating the values at risk 
  (u + (beta/xi) * ( ( ( (1 - alpha) / (Fbar) ) ^ (-xi) ) - 1) )
  
  }

Then we move on with the calculation of the value at risk for each observation in the test set. In order to do it, we will first create a vector to store the values, which will have a length of 300 as the number of observations in the test set. Then, we create a loop for the 300 observations, we want to create the window considered to calculate the value at risk by including the 365 observations that precedes the one for which we want to estimate the value at risk. Eventually, we get the different values at risk for each one of the observations.

#create a vector to store the values 
VAR <- as.numeric(1:300)

for(i in 1:300){
  #select window 
  wind <- sales %>% dplyr::slice((1913-i):(1913-i-365))
  #calculate VaR
  VAR[i] <- VaRest(wind)
}

VAR
##   [1] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
##   [9] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
##  [17] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
##  [25] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
##  [33] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
##  [41] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
##  [49] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
##  [57] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
##  [65] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6438
##  [73] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
##  [81] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
##  [89] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
##  [97] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
## [105] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
## [113] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
## [121] 138.6438 138.6504 135.7895 135.7767 135.7895 135.7895 135.7767 135.7767
## [129] 135.7767 135.7767 135.7767 135.7885 135.8013 135.8013 135.8013 135.8013
## [137] 135.8013 135.8013 135.8013 135.8013 135.8013 135.7886 135.7886 135.7886
## [145] 135.7886 135.7886 135.7886 135.7886 135.8015 135.8015 135.8015 135.8015
## [153] 135.8015 135.8015 135.8015 135.8015 135.8015 135.8015 135.8015 135.8015
## [161] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [169] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [177] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [185] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [193] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [201] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [209] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8310
## [217] 135.8310 135.8310 135.8310 135.8310 135.8310 135.8310 135.8310 135.8481
## [225] 135.8481 135.8481 135.8481 135.8481 135.8481 135.8481 135.8667 135.8667
## [233] 135.8667 135.8667 135.8667 135.8893 135.8893 135.8893 135.8893 135.8893
## [241] 135.8893 135.8893 135.8893 135.9120 135.9120 133.0342 121.1951 121.1566
## [249] 121.1566 121.1566 121.1566 121.1951 121.2387 121.5214 121.4619 121.5214
## [257] 121.5830 121.6408 150.9121 152.5583 152.4983 152.4983 152.4983 152.4983
## [265] 152.4953 152.4927 152.4927 152.4927 152.4927 152.4927 152.4927 152.4927
## [273] 152.4927 152.4904 152.4904 152.4904 152.4904 152.4904 152.4904 152.4904
## [281] 152.4904 152.4927 152.4927 152.4927 152.4927 152.4927 152.4904 152.4904
## [289] 152.4904 152.4904 152.4904 152.4904 152.4904 152.4904 152.4904 152.4927
## [297] 152.4927 152.4927 152.4927 152.4927

We can see that there are quite a lot of repetitions in the estimation of the values at risk, this is given by the fact that the the observations contained in the windows of adjacent dates are all the same but one, as we discard the oldest one to replace it with the real value of the observation of the day before the one we are considering. Hence, if none of these two values is describing a violation (meaning that it is higher than the threshold that we have chosen), the computation of the value at risk won’t be affected, as all the paramters that we use to estimate the VaR will remain unchanged: the weights will be the same (since we consider only the values that are above the threshold) and so will be the GPD and its parameters’ estimations, the fbar is also the same as Nu stays the same and the alpha is fixed.

3. Violations

Now, we compare the VaR we estimated to the real values, to see how many violations there are. The should be following a bernoulli distribution with n = 300 and p = 1 - alpha.

#number of observed violations
violations <- sum(ifelse(VAR > salesTest, 0, 1))
#number of theoretical violations, given by the mean of the bernoully distribution
ev <- 300*(1 - p)

We find out that there are 21 violations.

We can compare this value to the theoretical one that we expect to get, which is 27.2727273.

Binomial test

We can test whether they are statistically significantly different by using a binomial test, which has the following characteristics:

\[H_0: \pi = \pi_0\] In our case the theoretical π will equal to 0.0909091, while the probability found in the observations is 0.07. We want to test if waht it observed is different from the expected value.

(bintest <- binom.test(violations, 300, p = 1-p))
## 
##  Exact binomial test
## 
## data:  violations and 300
## number of successes = 21, number of trials = 300, p-value = 0.2287
## alternative hypothesis: true probability of success is not equal to 0.09090909
## 95 percent confidence interval:
##  0.04385034 0.10501400
## sample estimates:
## probability of success 
##                   0.07

We can see that the p-value is quite high (0.2286663), more specifically is higher than the significance level alpha = 5%, hence we cannot reject the null hypothesis of the two proportions being the same and hence that the distribution of the violations in the windows considered follows a binomial(300, 1-p).

In this third practical, we expand the retail sales dataset of Practical 2 to include multiple products. Data for Product 1 (from Practical 2) and an additional 9 products are given in two forms: \(sales10.csv\) and \(sales10adjusted.csv\), with the same format. As the name indicates, sales volumes in the latter file have been adjusted for trend and seasonality (removing growth and cyclical patterns).

#data

sales10 <- read.csv(here::here("sales_10.csv"))
sales10adj <- read.csv(here::here("sales_10_adjusted.csv"))

Exercise A

Explain the challenges involved when analysing the risk of stock-out across several products. What is the danger of considering each series independently?

The aplication of an independently analysis (EVT) could cause a loss of information and therefore a lack of preparation in face of a risky event. For example, the analysis of stock out from 3 products independently are not consider risky, nevertheless, the stock out of the 3 together could be consider as a very high risky event due to the high probability of loosing a customer.

MEVT has as a main objective to find the probabilities that extreme events happen together.

The main problem would be given by the fact that the different products may not be independent even if we consider them to be, it can be seen if one of them has an extreme value in the demand and so does another one, every time at the same date. This can happen as the products may be related, as for example could be bread and butter, meaning that if there is an increase in the demand of one product, so there will be in the demand for the other one. Another possibility is the fact that, especially for extreme values, so an abrupt change in the demand, could be determined by external factors, such as a natural disaster, that would make the demand for all sort of products increase from one day to the other. Indeed, it could be a problem, as by analysing the stock-outs, which are extreme values, independently one from the other, we could miss some important relationship that are between them and not being able to assess the risk in the best way possible.

Exercise B

Propose one graphical and one numerical method of detecting dependence of extreme values of the demand across several products. Apply your chosen methods to the adjusted sales data and identify groups of related products (if any exist).

n° Analysis Method
1 Graphical We will apply two methods:
1) A plot with all the products distributions over time, in which we will set a threshold as an arbitrary value, in order to enable us the study the ones which are above it and see if there is repeated pattern (for example having always two products that are over the threshold in two close dates)
2) A chi-plot among all the different couples of products, so that we can study the distribution of the tails of the products’ demands.
2 Numerical We will compute the chibar values, which will confirm us the tail dependences among all the different couples of products.
3 Aggregation In this step we will focus in the grouping of products, in order to do so, we will use the following 2 method:
1) We will calculate the mean and variance for each product and display them on a graph, to see if they have any other product close to them
2) Then, we will compute the Kendall-Taus distances among the products and see the clustering of the products based on it.

Graphical method

First, let’s create a plot with all the products’ demands distibutions.

All the products distribution

## integer(0)

In order to determine the threshold for the extreme values, we will take the 98% quantile of the distribution of the sales of each product, as we could find very different distributions and taking only one threshold overall would not give us any important information.

Chi-plot by couple of products

As we can see, we have some observations in each product that are above the line. We are particularly interested in seeing if there is a repetition in the observations that are above the lines that happen around the same time for different products.

Repitition Day

Now, we will determine the number of the observation in which we pass the line for each products and we try to see if there are some days that are close enough or repetead.

Table 2: Check if there are any repetition above the lines that happen around the same time between the products
Product 1 Product 2 Product 3 Product 4 Product 5 Product 6 Product 7 Product 8 Product 9 Product 10
240 300 205 66 299 8 209 242 102 39
299 301 301 154 330 9 223 253 104 53
329 329 329 301 630 325 325 321 226 101
336 332 332 332 663 330 349 334 249 193
337 382 382 333 696 337 365 345 282 226
664 387 387 336 940 522 380 372 283 248
696 665 388 584 946 664 413 413 307 255
703 696 576 665 947 696 468 436 333 273
938 888 664 694 976 703 471 611 437 282
960 933 696 696 977 945 497 618 499 283
962 1061 1061 699 1039 955 584 788 522 285
1018 1226 1281 888 1042 1038 701 974 603 317
1019 1274 1304 1060 1044 1042 758 1263 736 432
1039 1278 1426 1112 1063 1043 771 1264 744 470
1041 1279 1428 1252 1068 1044 779 1335 745 590
1044 1312 1526 1253 1085 1068 784 1339 829 614
1068 1399 1554 1275 1398 1136 828 1690 1166 645
1136 1426 1790 1399 1402 1191 918 1728 1261 683
1664 1428 1791 1425 1475 1526 1435 1736 1382 1166
1762 1791 1793 1428 1762 1763 1517 1910 1746 1382

We can see that especially for the dates of products 2 and 3 we have some repetitions of the same values, also for products 8 and 7 it looks like there is some similitudes or closeness of the values.

Chi-plot

The other method we can use is to create the chi plot for all the couples of products to assess the distributions of their tails and see if there is any correlation.

Key takeaway

We can see that in the highest quantiles we have a straight line for the majority of products, only for a couple of them (like p1p3 and p5p6), there seems to be a line that is not centered around zero, but to a higher value, meaning that there could be some correlation among them, which could be a problem.

Matematical method

While for the mathematical method, we will construct a matrix containing the chi bar values for all the couples of products.

(#tab:Kable of tail dependance)Chi bar Matrix
p2 p3 p4 p5 p6 p7 p8 p9 p10
p1 0.218 0.181 0.200 0.442 0.442 0.039 0.218 0.094 -0.072
p2 0.000 0.508 0.430 0.236 0.140 0.161 0.068 0.094 0.007
p3 0.000 0.000 0.406 0.218 0.218 0.118 0.094 0.007 0.068
p4 0.000 0.000 0.000 0.161 0.140 0.039 0.161 0.094 -0.072
p5 0.000 0.000 0.000 0.000 0.430 0.039 0.118 0.039 -0.029
p6 0.000 0.000 0.000 0.000 0.000 -0.029 -0.029 0.007 -0.127
p7 0.000 0.000 0.000 0.000 0.000 0.000 0.094 0.039 0.007
p8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.068 0.118
p9 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.328

We can see that the highest values appear to be for products 2 with product 3 and 4, product 1 with product 5 and 6, product 3 with 4, products 5 with 6. This is very much in line with we have found in the chi plot that we created above.

Aggregation method

To investigate a possible clustering of the different products, we will create create a graph with the mean and variance for each product.

Distance: Mean and variance

Defining k

Cluster with distance: Kendall’s tau

Key takeaway What is interesting to notice, is that the clustering are the same if we consider the whole distribution and only the most extreme values. Netherless, the correlation are stronger in the extreme values, thus we confirm that an independent analysis could potentialy lose important information.

In both cases, we can see that there is some similarity in some products, and some differences in others, especially product 1 appreas to be very different from the others, while products 9, 10 7, 8, 4, 3 and 2 appear very similar (the cluster in green). However, the number of clusters that should be used is not very clear, as the decrease of total within sum of sqaure is quite linear.

Exercise C

In view of your answer to (b), would you apply a Gaussian copula model to these data?

As well as the EVT the Gaussian model does not caputre all the information for the extreme values. Keeping this concept for multivariate extreme value theory we can see that in this case the corelation is higher, thus we cannot apply this model to fit the extreme value of the products’ distribution.

As we could find some depence in the distribution of the tails of our products’ demand, we believe that a gaussian distribution would not be optimal, as it is not able to consider the depences in the tails of the distribution.

Exercise D

Fit a multivariate model to the adjusted sales data and estimate the 95% Demand at Risk for the sum:

\[S = X^{(1)} + X^{(2)} + ... + X^{(10)}\]

  1. (Hard) Fit a multivariate model to the adjusted sales data and estimate the 95% Demand at Risk for the sum: \[S = X^{(1)} + X^{(2)} + ... + X^{(10)}\] You can do so with the following steps:
# These functions can be used to transform data to uniform scale,
# when the data are analysed separately below and above a threshold
# e.g. applying the GPD to values above a high threshold and treating
#      those below the threshold with the empirical CDF.

transform_to_uniform <- function(y,
                                 cdf,
                                 threshold,
                                 ...) {
  # ...: arguments to `cdf`
  
  
  above <- y > threshold
  # proportion of exceedances
  p <- mean(above)
  # under the threshold, use the empirical CDF
  # see ?ecdf
  ecdf_below <- ecdf(y[!above])
  empirical <- ecdf_below(y[!above])
  # above, apply cdf
  theoretical <- cdf(y[above], mean.sales10adj)
  
  transformed <- numeric(length = length(y))
  # "Glue together" the empirical and theoretical parts by rescaling them:
  
  #        empirical          e.g. N(0,1)
  # [-----------------------|-------------]
  # 0                      1-p            1
  
  transformed[!above] <- (1 - p) * empirical    # (1)
  transformed[above] <- 1 - p + p * theoretical # (2)
  
  return(list(
    transformed = transformed,
    ecdf = ecdf_below,
    prop = p
  ))
}

inverse_transform <- function(u, ecdf, quantile_function, p, ...) {
  # ...: arguments to quantile_function
  # Basically, undo the transformations (1-2) above
  above <- u > 1 - p
  original_scale <- numeric(length = length(u))
  original_scale[!above] <- quantile(ecdf, u[!above] / (1 - p))
  original_scale[above] <-
    quantile_function((u[above] - (1 - p)) / p, ...)
  return(original_scale)
}

  1. Transform the sales for each product to a uniform scale, that is, compute \[F_i(X^{(i)})\] for each X in 1, …, n being each of the product sales and F being the estimated Extreme Value distributions.

  2. Fit a multivariate copula of your choice using the copula::fitCopula function.

library(copula)
## Warning: package 'copula' was built under R version 3.6.2
unif <- cbind(unif1$transformed, unif2$transformed, unif3$transformed, unif4$transformed, unif5$transformed, unif6$transformed, unif7$transformed, unif8$transformed, unif9$transformed, unif10$transformed)

gum.fit <- gumbelCopula(5, dim=10)
clay.fit <- claytonCopula(4, dim=10)

fit.gum <- fitCopula(gum.fit, method = "itau", data=unif)
## Warning in .local(copula, tau, ...): For the Gumbel copula, tau must be >= 0.
## Replacing negative values by 0.
summary(fit.gum)
## Call: fitCopula(copula, data = data, method = "itau")
## Fit based on "inversion of Kendall's tau" and 1913 10-dimensional observations.
## Gumbel copula, dim. d = 10 
##       Estimate Std. Error
## alpha     1.04         NA

  1. Simulate from your fitted copula and transform the values back to their original scales (i.e. undo the transformation in i.)

  2. Compute the simulated values of S and their 95% Value at Risk. 

## 'data.frame':    50 obs. of  1 variable:
##  $ rev10: num  1.85 3.67 2.38 1.38 1.5 ...
rev1 rev2 rev3 rev4 rev5 rev6 rev7 rev8 rev9 rev10
1.833397 1.721064 1.674045 2.215008 1.289980 1.6627326 1.729129 1.754176 1.383327 1.848719
1.389896 2.381227 2.229203 1.342169 1.162100 1.1612029 2.654195 2.021832 1.939627 3.669373
1.613484 1.321973 1.447119 1.644108 1.619078 1.1096222 1.435127 1.564317 1.424517 2.384326
1.785663 1.704681 1.697891 1.674859 1.839798 2.1217804 1.311406 1.458367 1.444699 1.379314
1.208821 1.198224 1.581209 1.519660 1.101783 1.4266296 1.190012 1.490266 1.510652 1.496257
1.634306 1.590843 1.607376 1.675039 1.295762 1.7613302 3.214484 1.391533 1.703005 1.556532
1.533151 1.605448 1.601066 1.915172 1.079881 1.2236936 1.330645 1.609578 2.456184 1.364795
1.190154 2.004760 1.755000 1.918430 1.271633 1.2601710 1.193537 1.671078 1.876376 1.357561
1.233841 1.562571 1.978563 2.452361 1.281525 1.0810867 2.127564 1.349537 1.920688 1.536425
1.496956 1.982314 2.487770 1.444414 2.045698 1.2799271 2.037815 1.575224 1.249127 1.366893
2.113782 1.321333 1.474030 1.897923 1.423585 1.2549105 2.865743 1.721762 1.800719 1.202008
1.607884 2.760839 1.988598 1.933120 1.177907 2.4900720 1.756004 1.597384 1.611118 1.663281
2.873292 1.842880 1.519386 2.084405 1.643726 1.2696415 1.459083 1.387173 1.477006 1.200218
1.210557 1.101677 1.461695 1.356342 1.226587 2.0802226 1.816353 1.888337 1.265161 1.740980
1.585486 3.635026 1.416224 1.432849 1.169583 2.0454140 1.214606 1.546852 1.169102 1.778391
1.331183 1.123347 1.492611 2.248024 1.785575 1.7228191 1.311531 1.902817 1.256973 1.841823
1.256543 1.688945 1.376746 1.544800 1.525465 0.9615920 1.199900 1.537253 1.552004 1.660870
1.296168 1.527678 1.485815 1.557216 1.664586 1.9478886 2.309075 1.417349 1.952647 1.339120
1.181620 2.049961 2.321236 1.355867 1.337488 2.2613678 1.656007 1.377070 1.485760 1.945392
2.100210 1.077789 1.433739 1.592683 1.418845 1.9574048 1.806083 1.486208 2.004329 1.487589
2.061404 1.913821 1.526011 1.564677 1.936295 1.6178485 1.749009 2.314871 2.745721 1.265690
1.517808 2.310673 1.970856 1.690513 2.531249 1.1877103 1.629864 1.716186 1.393269 1.860563
1.733879 2.434929 2.002785 1.615377 1.387149 1.9369006 2.445183 1.914617 1.360017 2.552996
1.089229 1.686079 1.450542 2.046041 2.275769 1.7553744 1.859622 1.632278 1.644185 1.424989
1.442043 1.457383 1.398268 1.836652 1.536172 0.9777846 1.261231 1.851284 1.188082 2.038215
1.659330 1.776774 1.362589 1.640176 1.108892 1.2932010 1.259842 2.013255 1.482082 1.527087
2.070678 1.731888 2.068902 1.868487 1.719934 1.8460408 1.547536 2.343498 2.364856 1.509815
1.474110 1.288580 1.660731 1.388430 2.238346 1.0857424 1.239393 1.507835 1.746296 1.285536
1.352808 1.179803 2.185961 1.485179 1.554753 1.1317846 2.322521 1.376496 2.409522 1.659858
2.278935 1.112572 1.419617 1.421075 1.114297 1.6818651 1.509684 1.583086 1.365620 1.382108
1.370466 2.033608 1.692362 1.491893 1.742530 1.4601652 1.220284 1.840711 1.532514 1.341479
2.271847 1.241626 2.203395 1.500229 1.344072 0.9984627 1.136363 1.975560 1.423543 1.359376
1.275806 1.911564 2.329551 1.459843 1.085427 1.3701531 1.320556 2.467837 2.056330 1.359570
2.040185 1.360157 1.344735 2.155631 1.550392 1.1173153 1.438916 1.965333 1.369547 2.068387
1.358254 1.464846 2.032040 1.799191 1.976331 1.2252450 1.448942 1.272767 1.345834 2.668137
1.214027 1.105772 1.565685 2.043336 1.645061 1.8596732 1.326358 1.691640 1.362605 2.307718
1.831860 2.649468 1.902648 1.695757 2.128824 2.8308023 1.248845 1.895207 1.390897 1.223339
1.946117 1.639319 1.518642 1.798334 1.470069 2.4631674 1.379029 2.494113 1.425435 1.330472
1.647007 1.389182 1.393164 1.372668 1.168176 1.4802188 1.445032 1.802889 1.339266 1.513919
2.111566 1.621631 2.088837 1.377249 2.392124 1.0427114 3.191552 1.875108 2.859082 1.236406
2.000295 1.516690 1.353543 2.443163 1.752450 3.0232093 1.307114 1.446175 1.901184 1.478862
1.430779 1.093592 1.640822 1.374975 1.842401 1.7282744 1.288309 1.915499 1.637622 1.528434
2.191231 1.139949 2.386237 2.314526 1.615161 1.1275522 1.909695 2.169519 1.483068 1.271424
2.358958 1.933521 2.134471 1.374123 1.654963 2.2836139 1.991499 1.647774 2.126623 1.221293
1.525797 2.363068 1.526451 1.380267 2.215259 1.9524406 1.788769 1.877078 1.670259 1.261554
1.901465 1.395117 1.653835 1.401864 1.668551 2.5688082 1.586492 1.833891 2.193103 1.581098
1.454989 1.165118 1.545595 2.484355 2.508851 1.2133469 1.145841 1.721651 1.805971 1.391955
1.721983 1.141894 3.133116 1.607453 1.309691 0.9835133 1.310715 2.029665 1.703807 2.299050
2.038518 1.195019 2.221920 1.631545 1.242957 1.3963297 1.266906 1.292525 1.666654 1.737291
1.322076 1.525024 1.427202 1.531989 1.465071 1.4171284 1.581289 1.596690 1.173473 1.476023
Table 3: Extreme Values Data frame
rev1 rev2 rev3 rev4 rev5 rev6 rev7 rev8 rev9 rev10
1.833397 1.721064 1.674045 2.215008 1.289980 1.6627326 1.729129 1.754176 1.383327 1.848719
1.389896 2.381227 2.229203 1.342169 1.162100 1.1612029 2.654195 2.021832 1.939627 3.669373
1.613484 1.321973 1.447119 1.644108 1.619078 1.1096222 1.435127 1.564317 1.424517 2.384326
1.785663 1.704681 1.697891 1.674859 1.839798 2.1217804 1.311406 1.458367 1.444699 1.379314
1.208821 1.198224 1.581209 1.519660 1.101783 1.4266296 1.190012 1.490266 1.510652 1.496257
1.634306 1.590843 1.607376 1.675039 1.295762 1.7613302 3.214484 1.391533 1.703005 1.556532
1.533151 1.605448 1.601066 1.915172 1.079881 1.2236936 1.330645 1.609578 2.456184 1.364795
1.190154 2.004760 1.755000 1.918430 1.271633 1.2601710 1.193537 1.671078 1.876376 1.357561
1.233841 1.562571 1.978563 2.452361 1.281525 1.0810867 2.127564 1.349537 1.920688 1.536425
1.496956 1.982314 2.487770 1.444414 2.045698 1.2799271 2.037815 1.575224 1.249127 1.366893
2.113782 1.321333 1.474030 1.897923 1.423585 1.2549105 2.865743 1.721762 1.800719 1.202008
1.607884 2.760839 1.988598 1.933120 1.177907 2.4900720 1.756004 1.597384 1.611118 1.663281
2.873292 1.842880 1.519386 2.084405 1.643726 1.2696415 1.459083 1.387173 1.477006 1.200218
1.210557 1.101677 1.461695 1.356342 1.226587 2.0802226 1.816353 1.888337 1.265161 1.740980
1.585486 3.635026 1.416224 1.432849 1.169583 2.0454140 1.214606 1.546852 1.169102 1.778391
1.331183 1.123347 1.492611 2.248024 1.785575 1.7228191 1.311531 1.902817 1.256973 1.841823
1.256543 1.688945 1.376746 1.544800 1.525465 0.9615920 1.199900 1.537253 1.552004 1.660870
1.296168 1.527678 1.485815 1.557216 1.664586 1.9478886 2.309075 1.417349 1.952647 1.339120
1.181620 2.049961 2.321236 1.355867 1.337488 2.2613678 1.656007 1.377070 1.485760 1.945392
2.100210 1.077789 1.433739 1.592683 1.418845 1.9574048 1.806083 1.486208 2.004329 1.487589
2.061404 1.913821 1.526011 1.564677 1.936295 1.6178485 1.749009 2.314871 2.745721 1.265690
1.517808 2.310673 1.970856 1.690513 2.531249 1.1877103 1.629864 1.716186 1.393269 1.860563
1.733879 2.434929 2.002785 1.615377 1.387149 1.9369006 2.445183 1.914617 1.360017 2.552996
1.089229 1.686079 1.450542 2.046041 2.275769 1.7553744 1.859622 1.632278 1.644185 1.424989
1.442043 1.457383 1.398268 1.836652 1.536172 0.9777846 1.261231 1.851284 1.188082 2.038215
1.659330 1.776774 1.362589 1.640176 1.108892 1.2932010 1.259842 2.013255 1.482082 1.527087
2.070678 1.731888 2.068902 1.868487 1.719934 1.8460408 1.547536 2.343498 2.364856 1.509815
1.474110 1.288580 1.660731 1.388430 2.238346 1.0857424 1.239393 1.507835 1.746296 1.285536
1.352808 1.179803 2.185961 1.485179 1.554753 1.1317846 2.322521 1.376496 2.409522 1.659858
2.278935 1.112572 1.419617 1.421075 1.114297 1.6818651 1.509684 1.583086 1.365620 1.382108
1.370466 2.033608 1.692362 1.491893 1.742530 1.4601652 1.220284 1.840711 1.532514 1.341479
2.271847 1.241626 2.203395 1.500229 1.344072 0.9984627 1.136363 1.975560 1.423543 1.359376
1.275806 1.911564 2.329551 1.459843 1.085427 1.3701531 1.320556 2.467837 2.056330 1.359570
2.040185 1.360157 1.344735 2.155631 1.550392 1.1173153 1.438916 1.965333 1.369547 2.068387
1.358254 1.464846 2.032040 1.799191 1.976331 1.2252450 1.448942 1.272767 1.345834 2.668137
1.214027 1.105772 1.565685 2.043336 1.645061 1.8596732 1.326358 1.691640 1.362605 2.307718
1.831860 2.649468 1.902648 1.695757 2.128824 2.8308023 1.248845 1.895207 1.390897 1.223339
1.946117 1.639319 1.518642 1.798334 1.470069 2.4631674 1.379029 2.494113 1.425435 1.330472
1.647007 1.389182 1.393164 1.372668 1.168176 1.4802188 1.445032 1.802889 1.339266 1.513919
2.111566 1.621631 2.088837 1.377249 2.392124 1.0427114 3.191552 1.875108 2.859082 1.236406
2.000295 1.516690 1.353543 2.443163 1.752450 3.0232093 1.307114 1.446175 1.901184 1.478862
1.430779 1.093592 1.640822 1.374975 1.842401 1.7282744 1.288309 1.915499 1.637622 1.528434
2.191231 1.139949 2.386237 2.314526 1.615161 1.1275522 1.909695 2.169519 1.483068 1.271424
2.358958 1.933521 2.134471 1.374123 1.654963 2.2836139 1.991499 1.647774 2.126623 1.221293
1.525797 2.363068 1.526451 1.380267 2.215259 1.9524406 1.788769 1.877078 1.670259 1.261554
1.901465 1.395117 1.653835 1.401864 1.668551 2.5688082 1.586492 1.833891 2.193103 1.581098
1.454989 1.165118 1.545595 2.484355 2.508851 1.2133469 1.145841 1.721651 1.805971 1.391955
1.721983 1.141894 3.133116 1.607453 1.309691 0.9835133 1.310715 2.029665 1.703807 2.299050
2.038518 1.195019 2.221920 1.631545 1.242957 1.3963297 1.266906 1.292525 1.666654 1.737291
1.322076 1.525024 1.427202 1.531989 1.465071 1.4171284 1.581289 1.596690 1.173473 1.476023